Comparison of colostrum and milk extracellular vesicles small RNA cargo in water buffalo

Recently, much interest has been raised for the characterization of signaling molecules carried by extracellular vesicles (EVs), which are particularly enriched in milk (mEVs). Such interest is linked to the capability of EVs to cross biological barriers, resist acidification in the gastric environment, and exert modulation of the immune system, mainly through their microRNA (miRNA) content. We characterized the small-RNA cargo of colostrum EVs (colosEVs) and mEVs from Italian Mediterranean buffalo through next generation sequencing. Colostrum (first milking after birth) and milk (day 50 of lactation) were sampled from seven subjects from five farms. ColosEVs and mEVs were subjected to morphological characterization, followed by high-depth sequencing of small RNA libraries produced from total RNA. The main difference was the amount of EV in the two samples, with colostrum showing 10 to 100-fold higher content than milk. For both matrices, miRNA was the most abundant RNA species (95% for colosEVs and 96% for mEVs) and three lists were identified: colosEV-specific, mEV-specific and shared most expressed. Gene ontology (GO) enrichment analysis on miRNA targets highlighted many terms related to the epigenetic, transcriptional and translational regulations across the three lists, with a higher number of enriched terms for colosEV-specific miRNAs. Terms specific to colosEVs were related to “cell differentiation” and “microvillus assembly”, while for mEV “cardiac and blood vessel development” and “mitochondria” emergerd. Immune modulation terms were found for both sample-specific miRNAs. Overall, both matrices carry a similar molecular message in terms of biological processes potentially modulated into receiving cells, but there is significant difference in the abundance, with colostrum containing much more EVs than milk. Moreover, colosEVs carry molecules involved in signal transduction, cell cycle and immune response, as for mEVs and EVs of other previously characterized species, but with a special enrichment for miRNAs with epigenetic regulation capacities. These beneficial characteristics of colosEVs and mEVs are essential for the calf and could also be exploited for the therapeutic purposes in humans, although further studies are necessary to measure the sanitization treatment impact on EV conservation, especially in buffalo where milk is consumed almost exclusively after processing.

Estrogen receptor α RUNX1T1 Runt-related transcription factor 1, translocated to 1 MED Mediator complex DRIP Vitamin D receptor-interacting protein TIF1A Transcription intermediary factor 1 α Milk constitutes the main and mostly unique source of nutritional elements for mammalian newborns, and even more, a complex delivery system for maternal messages toward the offspring 1 .The presence of bioactive compounds in milk makes it an important source of nutrients for human consumers as well 2 .Milk is richer in all major constituents, with high energy and nutritional value conferred by fat, which is one of the main fractions 3,4 .Moreover, buffalo milk offers advantages over milk from other species in promoting health due to specific properties such as a higher protein content, and calcium, the main mineral, 1.5-fold higher than in cow milk 3 .Furthermore, milk from buffalo is characterized by a high content of antioxidant and anti-inflammatory molecules such as tocopherols and vitamin A 2,4 , and contains exclusive compounds including biliverdin, bioactive pentasaccharides, and gangliosides 5 .The microbial composition of buffalo milk has also been evaluated in detail, identifying characteristics that qualify it as a functional food with probiotic effects 6 .As expected, the milk composition of these bioactive molecules can change over time, in relation to the lactation period, environmental stimuli, diet, and genetic factors [7][8][9][10] .
Colostrum shows a significant different composition compared to milk 11 ; for buffaloes, as for all other ruminant species, colostrum represents a crucial source of immunoglobulins for newborns since the syndesmochorial placenta inhibits the transfer of antibodies.This aspect makes buffalo newborns particularly susceptible to infections, and it is the main reason for the importance of colostrum intake 12 .In addition to the higher protein content compared with mature milk, colostrum also contains greater levels of total solids, non-fat solids, and milk urea nitrogen, while exhibiting lower levels of fat and lactose 13 .For minerals, buffalo colostrum has high concentrations of Na, Mg, Co, Fe, and K, while milk is more concentrated in Ca 13 .
Despite the importance of buffalo milk for dairy products, few studies have been conducted to evaluate its colostrum composition, and only one evaluated buffalo milk extracellular vesicles (mEVs) 14 .Extracellular vesicles (EVs) are spherical micro/nano-sized cell structures 15,16 characterized by a phospholipid bilayer enclosing a plethora of molecules, including RNAs, lipids, proteins and DNA fragments [17][18][19][20] .Basing on recent developed guidelines of the international society for extracellular vesicles (ISEV) 21 and the increasing knowledge about these particular structures, EVs are mainly subdivided into three categories with respect to the size: small (30-150 nm), medium-large (100-1000 nm) and apoptotic bodies (50-5000 nm).Small and medium-large EVs can refer to two main vesicle typologies, exosome and microvesicles, which differ in the biogenesis process 22 .Since no specific molecular markers exist to distinguish between these two subtypes and the co-isolation quite always occurs, it is desirable to use the nomenclature based on size (small/medium-large EVs).
EVs can be released by all cell types in the extracellular matrix and can reach close and distant receiving cells, which, after the vesicle absorption, can undergo the modulation of cellular processes 23 .Indeed, EV uptake can induce a modification in the recipient cell toward an anti-inflammatory, anti-cancer, pro-regenerative or angiogenetic phenotype.These modifications depend on various factors including the recipient cell type, the source cell and the environmental stimuli received 15,24,25 .These effects are mediated by the molecular cargo, which is protected within EVs, thus keeping the integrity and functionality, particularly concerning microRNAs (miRNAs) 20 .The biocompatibility and the protection conferred by EVs to the cargo, allow these molecules to overcome biological barriers and adverse environments, such as the gastrointestinal tract, without damage, thus maintaining their functionality 26 .A study conducted in mice, reported that miRNAs loaded in mEV, can be detected in distant tissues after vesicle ingestion 27 .

Colostrum and milk collection and extracellular vesicle isolation
From each of seven Italian Mediterranean buffalos, reared in five different farms located in central Italy (Lazio region), 500 ml from the whole milking of colostrum (first milking) and milk (fifty days after parturition) were obtained.Rearing conditions are standard intensive livestock systems.The seven subjects calved within 4 days range (from January 24th to January 28th, 2022).Individual colostrum and milk samples, from the same seven buffalos, were packaged in 50 ml plastic tubes containing Bronopol ® (2-bromo2nitropropano1,3diol) as a preservative, and immediately stored at 4 °C and processed within 24 h, avoiding cryo-preservation to minimize artifacts.
To isolate colostrum EVs (colosEVs) and mEVs, the protocol reported by Mecocci et al. 38 was used.Colostrum and milk EVs were isolated through serial differential centrifugations (DC), and serum was treated with ethylenediaminetetraacetic acid tetrasodium salt dihydrate (EDTA) to precipitate most of the contaminating proteins.Fat globules and cellular debris and protein complexes were removed throw by applying preliminary DC steps: two sequential 3000 × g centrifugations for 10 min at room temperature (Eppendorf ® Centrifuge 5810R with a F34-6-38 rotor) and then a 10,000 × g for 1 h at 4 °C after the incubation for 15 min in ice with an equal volume of 0.25 M EDTA (pH 7.4).Afterward, ultracentrifugation at 35,000 × g for 1 h at 4 °C (Optima L-100 XP, Beckman Coulter, Milano, Italy) with a Type 45 Ti rotor (Beckman coulter) was carried out to eliminate precipitated proteins from sera, and a final ultracentrifugation at 200,000 × g for 90 min at 4 °C was applied to recover EV pellets.

Colostrum and milk EV characterization
The effective isolation of colosEVs and mEVs was assessed by Transmission Electron Microscopy (TEM) observation and ExoView™ R100 technology (NanoView Biosciences, Brighton, MA, USA).
A drop of mEV suspension from one pellet was placed on Parafilm.A Formvar-coated copper grid (Electron Microscopy Sciences) was placed over each drop for approximately 30 min to allow the mEVs to adhere to the surface.The grids were then washed in PBS and distilled water and then contrasted with 2% uranyl acetate for 5 min.The samples were observed using a Philips EM208 transmission electron microscope equipped with a digital camera (University Centre of Electron and Fluorescence Microscopy-CUMEF).The ExoView™ R100 (NanoView Biosciences, Brighton, MA, USA) allowed us to test the isolated EVs for the positivity to EV markers such as cluster of differentiation 81 (CD81), CD9 and CD63, on a chip through an antibody-mediated capture.First, colosEV and mEV pellets were resuspended in 400 μl of 1 × PBS, filtered (0.22 μm pore size) and diluted 50 × in the incubation solution buffer of the human tetraspanin cargo kit (NanoView Biosciences, Brighton, MA, USA).Then, 50 μl of each solution were allowed to capture on the chip, left in incubation and washed, and finally made to react against the same fluorescently marked antibodies (CD81, CD9 and CD63).This technic also allowed to measure the size distribution and concentration of EV solutions.

RNA extraction and library preparation
From each sample, four EV pellets were treated with 1 ml (each) of TRIzol™ (Thermo Fisher Scientific, Waltham, MA, USA) in order to extract total RNA content using the miRNeasy Mini Kit (QIAGEN, Germantown, MD, USA) and following the manufacturer's instructions.On-column DNase digestion (RNase-Free DNase Set, QIAGEN, Germantown, MD, USA) was applied and the each extracted RNA (from 7 colosEV and 7 mEV Vol:.( 1234567890

Bioinformatic analysis
Raw reads from Illumina sequencer were checked for quality through the FastQC tool (https:// www.bioin forma tics.babra ham.ac.uk/ proje cts/ fastqc/) and low quality/short reads and adaptors trimmed using TrimGalore (https:// github.com/ Felix Krueg er/ TrimG alore).BowTie2 51 was used to align cleaned reads against reference features, setting parameters for short sequences.In particular, a two-step mapping procedure was carried out, first aligning on Bos taurus miRBase database (v.22) hairpin micro RNA (miRNA) sequences 52 and then, the unmapped reads were used on Bos taurus reference genome (ARS-UCD1.2) to retrieve accurate information on miRNAs and other small RNAs.Unfortunately, we could not use the Bubalus bubalis genome and the specific miRBase since this species lacks annotation on miRNAs, which are the most represented features in the EV cargo.The uniquely mapped reads where selected and a comprehensive count matrix was built after the read counting for each feature type through FeatureCounts 53 software, merging the information generated from the two alignments and summing counts for the same miRNA if covered by reads from both mapping steps.The normalized count matrix was generated through the DESeq2 package 54 in R environment (v.4.3.0) 55and used to evaluate feature abundances considering all genes with normalized counts greater than 1 in 4 out of 7 samples, for both colosEVs and mEVs.
An exploratory analysis was carried out on the same normalized count matrix by applying a principal component analysis on the top 100 features (plotPCA function of DESeq2), and a heatmap (pheatmap package).Finally, the differential gene expression analysis was carried out on genes with normalized counts greater than 1 in at least 7 out of 14 samples and considering as differentially expressed genes (DEGs) those with a False Discovery Rate (FDR) < 0.05 and a log 2 Fold Change (FC) less than − 1 or greater than 1 (|log 2 FC|> 1) in colosEVs compared to mEVs.The R package ggplot2 56 was used to produce a volcano plot.
To better evaluate a possible sample specific message, miRNAs with a good expression level (comprising the 99.9% of total normalized counts) were kept, eliminating scarcely represented features, and specific (i.e.present in only one kind of sample) miRNAs of colosEVs (colosEV-specific miRNAs) and mEVs (mEV-specific miRNAs) were identified drawing a Venn diagram (https:// bioin forma tics.psb.ugent.be/ webto ols/ Venn/).From this analysis a large core with miRNAs shared between the two conditions was highlighted and further refined by filtering for expression levels (top miRNAs comprising 95% of total normalized counts), in order to identify a set of miRNAs particularly enriched and shared between colosEVs and mEVs, ideally representing a buffalo specific signature (third column of Table 3).

Functional analysis
The functional analysis was carried starting from the three lists of miRNAs already described colosEV-specific; mEV-specific and core miRNAs (shared by colosEV and mEV groups).Target genes were identified utilizing MiRWalk 3.0 (http:// mirwa lk.umm.uni-heide lberg.de/), with the input being the name of the most sequenced miRNA (5p or 3p) obtained from literature via miRBase (v.22).As a result, a list of target genes for each miRNA relative to the binding site on the targeted mRNA (3′-UTR, 5′-UTR, or CDS-coding region) was obtained.Then, for a more robust functional analysis, the targets were filtered for the number of miRNA hits, focusing on those genes targeted by about 50% of input miRNAs.This filter aims to identify the most probable targets likely to be modulated in receiving cells by colosEV or mEV upon uptake.Finally, the three lists of target genes derived from the initial groups of miRNAs were analysed with the Cytoscape 3.9.1 suite 57 .A Protein-Protein Interaction Network (PPI) was constructed based on the IMEx database 58 , and the network was clustered according to the number and type of connections between the nodes using the clusterMaker 2.0 app 59 with the gLay option.With the most interconnected clusters (with 50 minimum number of interactions), a Gene Ontology (GO) enrichment analysis for biological process vocabulary was conducted, using the ClueGO application 60 .Results included only GO terms with a FDR < 0.05 (Benjamini Hockberg correction 61 ).Moreover, to ease the reader's experience, a summary for the enriched GO terms for each cluster was provided, indicating functional group names within each cluster.Functional groups are identified by the ClueGO software based on related and interconnected enriched GO terms for a particular biological function and named with the most enriched term (lowest p-value). of the scope of Directive 2010/63/EU (art.1.5.f"practices not likely to cause pain, suffering, distress or lasting harm equivalent to, or higher than, that caused by the introduction of a needle in accordance with good veterinary practice").

Extracellular vesicle evaluation
The presence of EVs has been demonstrated through the Exoview technology.We observe a difference in terms of quantity and dimensions for milk and colostrum preparations.Indeed, from the anti-CD9 capture spot, colosEVs were about 4.5 times more concentrated with a slightly bigger mean diameter (Fig. 1A).All samples showed a higher cross-reactivity with the anti-CD9 antibody compared to anti-CD81 reaction, both detected by interferometric imaging.Differently, in fluorescence, the signal was specific and significant only for the anti-CD9 antibody.It was impossible to obtain data relating to anti-CD63 spots, as the results obtained were below the minimum detection limit, as expected, since no cross-reactions against this antibody were shown for mEVs in our previous characterizations 35,38 .TEM analysis (Fig. 1B) revealed no specific differences between colosEVs and mEVs.The observed corpuscular elements were not particularly well preserved, but their morphology was consistent with typical EVs.They were between 40 and 300 nm in size, with most particles having a diameter of 100 nm or less.The field appeared to be fairly clean, apart from the presence of some filamentous structures, particularly in the colostrum.

Sequencing output
Sequencing results for the seven colosEVs (C) and mEV (M) samples are reported in Tables 1 and 2 for first and second mapping, respectively.In detail, a deep sequencing was carried out with over 33 million raw reads for sample on average which become about 28 million after trimming procedures and quality check.After the first alignment step on miRNA database, about 16% of the cleaned reads were uniquely aligned (Table 1), in line with previous experiments on mEVs 38 .
About 1.6 million of the unmapped reads were identified in the second alignment step on the Bos taurus genome.Comparable results were obtained for the two sample types (Table 2).

Small RNA cargo
As expected from the alignment results and from previous studies 38,64 , most of the reads aligned to the miRNAs, included more than 95% of normalized counts for colosEVs and 96% for mEVs (Fig. 2).A principal component analysis carried out on the count matrix shows a clear separation of the two groups along the first eigenvector (PC1), which explains 86% of the variance (Fig. 3A).This is confirmed by the heatmap analysis, which highlights two clearly defined clusters separated by sample type rather than individual animals (Fig. 3B).
This finding was further confirmed by the consistent number of differentially expressed features (Additional file 1).In detail, a total of 1504 differentially expressed genes (DEGs, log2 Fold Change − log2FC >|1| and adjusted p < 0.05) was found, 961 up-regulated and 543 down-regulated in colosEVs compared to mEVs (Fig. 4A; Additional file 1).Interestingly, an opposite scenario emerged for the predominant class miRNA: out of 212 differentially expressed features fewer (28) were up-regulated whereas a higher number of down-regulated (184) miRNAs in colosEVs compared to mEVs was found.Moreover, from the total miRNAs found in colosEVs (145) and mEVs (179), a consistent number was found expressed in a specific sample type (10 colosEV-specific and 44 mEV-specific) and completely absent (or with negligible counts) in the other type, as visible in the Venn diagram (Fig. 4B-first and second columns of Table 3).
These were considered as sample type specific miRNAs and used for the functional analysis, in order to identify colostrum and milk specific messages.A great number of miRNA, however, was present in both EV types (135, Fig. 4B).From these, a restricted list comprising the most expressed miRNAs in both sample types was generated (Table 3, third column) and used for the functional analysis, to highlight main shared functions for colosEVs and mEVs.Mean and standard deviation values of normalized counts for miRNAs of each of the three lists are reported in Additional file 2.

Functional analysis
The functional analysis was carried out as follows: targets were retrieved for each list of miRNAs (specific colosEVs, mEVs, and shared miRNAs), generating three lists of targets; for each one of the three lists of targets, a PPI-network was built and reorganized in clusters, identifying hub genes (central nodes); for each cluster, a functional analysis for enriched biological processes was carried out.In detail, for each miRNA list reported in Table 3, all the targets were retrieved and filtered on the number of miRNA hits, considering only those bound by 50% of input miRNAs (Table 4).Complete information for each target about binding miRNAs and binding site are reported in Additional file 3.
The functional analysis was carried out starting from filtered targets, generating a Protein-Protein Interactions (PPI) network through the IMEx database.From the 28 targets of specific colosEV miRNAs, a complex network with 3275 nodes and 3918 edges was built, while for specific mEV miRNAs 2139 nodes and 2385 edges were found from the 35 targets.For the most expressed shared miRNAs, a network with 2483 nodes and 2564 edges was obtained from the 18 targets.These networks were separately organized in clusters based on the number of interactions between proteins, as reported in Fig. 5, resulting in 13, 10 and 11 clusters for colosEV-specific (Fig. 5A), mEV-specific (Fig. 5B) and core miRNA (Fig. 5C) targets, respectively.
A total of 24, 14 and 13 hub protein-coding genes, corresponding to proteins with the highest number of interactions with other proteins (over 50) sharing similar biological functions (red dots, Fig. 5), were found for colosEV, mEV and core miRNA targets, respectively.Two transcription factors, i.e.SP1 and TF1A, are among them but others are related to transcription such as the transcriptional cofactors ARC70, the coactivator SNW1, the coregulator BCLAF3 or the mediators ARC77, MED4, MED12, MED30, and RGR1.It is worth noting that some node proteins (ARC34, ARC70, ARC77, KIAA0928, RGR1, SP1 and SUR2) are shared between targets of specific colosEV and mEV miRNAs (Fig. 5A,B), while APOB is a central node of both mEV miRNA and core mRNA targets (Fig. 5B,C).
For each of these clusters, a GO enrichment analysis for the "biological process" vocabulary was carried out through the ClueGO app and complete results are reported in Additional files 4-6 (if no data is reported, it means no enriched GO terms of functional group emerged).As a summary, Tables 5, 6 and 7 report the first 15 (lower adjusted p-value) functional groups in which are enclosed the enriched biological processes for each cluster.Functional groups are identified by the ClueGO software based on related and interconnected terms for a particular biological function among the enriched and are named as the most enriched comprise term (lower p-value).

Discussion
The findings of the present study on the small RNA content of buffalo colostrum and milk EVs might be interesting both for the development of the immune system and calf welfare.Moreover, they could shed light on the potential effect of buffalo milk could have on human nutrition.It is well known, indeed, that breast colostrum and milk have a key role in immune system enhancement and systemic disease resistance in infants; furthermore, breastfeeding is associated with reduction of morbidity and mortality 65 .Recently, strong immune regulatory functions have been recognized for the abundant immunoregulatory miRNAs in colosEVs and mEVs 64 .In particular, miRNAs are characterized by a strong evolutionary conservation across different animal species, leading to the possibility of interspecies cross-talk which potentially induces the regulation of several cellular processes 31,[48][49][50] .
The first interesting result of the EV characterization was the one obtained through the ExoView™Assay, which represents an implementation of the Nanoparticle Tracking Analysis (NTA) to evaluate EV concentration and size by using specific EV antigens for vesicle detection.This method has an important advantage over NTA, allowing the measurement of only EVs and not of all the suspended nanoparticles.On the other hand, it can be challenging to find the right cross-reactivity against species-specific antigens, particularly in lesser-studied species such as buffalo.Nevertheless, it has been possible to have considerable EV detection with the anti-CD9 antibody.CD81, www.nature.com/scientificreports/differently, showed a lower efficiency of capture and CD63 gave unreliable results, as already reported for cow mEVs 35,38 .As shown in Fig. 1, EVs were detected in both milk and colostrum preparations, although showing some differences in terms of quantity and dimension.As we expected, the colostrum was more concentrated in EVs than milk and was characterized by slightly larger vesicles.It is known that colostrum has considerable influence and importance for the newborn and is qualitatively different from milk.Indeed, colostrum not only possesses a greater quantity of proteins, mainly antibodies 12 , compared to milk, but also a different composition in EVs carrying their immunomodulatory message, as highlighted in this study.
Concerning the molecular cargo, the sequencing produced over 33 M raw sequences (reads) on average, a notably high number for small RNA investigations, which facilitate the detection of the great part of RNAs enclosed in colosEVs and mEVs.As expected from our previous study on milk-derived EV cargo characterization 38 , about 21.5% of quality-controlled reads are uniquely aligned to miRBase and genome, overall.Moreover, the paired sampling lowered the individual variability increasing the statistical power of data obtained.This allowed us to discern a distinct molecular cargo composition between colosEVs and mEVs as illustrated by the principal component analysis (PCA) and the heatmap (Fig. 3), which reflects the differences of the two matrices.Indeed, the principal component 1 (PC1) of panel A completely divides the two types of EVs, explaining 86% of total variability, and panel B shows two clear clusters deriving from the two sample matrices.It is worth noting that in this study, we obtained samples from the same seven animals in the two times.Hence, the PCA analysis of Fig. 3 also demonstrates that individual animals have minimal influence on the two profiles.
Most of the features carried by both colosEV and mEV cargos (about 95% of the total feature types) are represented by miRNAs and, due to their predominance and impact on the post-transcriptional gene expression regulation, we focused our attention on this class of RNAs, retrieving the most probable interested target genes and discussing related biological processes that can be affected in receiving cells once EVs are taken up.With an opposite trend compared to total DEGs, miRNAs showed to be more down-regulated (184) and less up-regulated (28) in colosEVs vs mEVs, identifying a sort of enrichment for miRNAs in mEV cargo.Moreover, we identified sample specific miRNAs, expressed in a sample type and completely absent (or with negligible counts) in the other (Venn diagram of Fig. 4B, Table 3).However, most of the miRNAs enclosed in colosEVs and mEVs are shared, in line with literature results on human milk 66 .In general, it is possible to identify some miRNA families particularly represented in both milk and colostrum EVs.Indeed, as for other species, btamir-30 and bta-mir-200 family members were found among the most abundant and shared between colosEVs and mEVs (Table 3) while a completely different result was obtained for other known miRNA families 29,67,68 .For instance, bta-mir-148 members were found in our buffalo colosEVs and mEVs although at concentrations under the threshold we applied to find the shared most expressed miRNAs, while for the let-7 family, generally hugely abundant in mEVs 38,69 , we found only bta-let-7f in both EVs at quite low concentrations.
Despite the apparent dilution of miRNAs within colosEV cargo, the functional analysis on their targets is characterized by a network with a much greater number of interactions than mEVs.Also, the number of hub genes from clusters and the enriched biological processes found statistically significant were greater for colosEVs than for mEVs, apparently indicating a higher degree of diverse actions associated with colostrum EVs (compare Tables 5 and 6, Additional file 4 and Additional file 5).
Interestingly, many GO terms related to the epigenetic regulation were found enriched in all the three lists of genes, although colosEV miRNA targets showed a higher number of enriched biological processes related to DNA metabolism and chromatin organization.Similar results were obtained for terms related to transcriptional and translational regulations including protein phosphorylation, being colosEV miRNA target list particularly interested by an enrichment of these processes.
Enriched GO terms related to cell cycle, cytoskeleton organization, vesicular transport and ion transport were found for targets of core shared miRNAs (Table 7 and Additional file 6) although further regarding colosEV miRNA targets.Interestingly, biological processes of neuronal development as well as protein kinase activity were found for all the three lists with tight-junction development found among core miRNAs.However, specific terms associated with cell differentiation and microvillus assembly were observed for colosEV miRNA targets.Concerning mEV miRNA targets, biological processes for cardiac and blood vessel development as well as mitochondria emerged.
Many biological processes related to the immune modulation were found for both colosEV and mEV samples: regulation of type I interferon and interleukin 1 (IL-1) production.Moreover, colosEV additionally showed "regulation of cytokine signalling pathway", "regulation of NFKB signal", "macrophage chemotaxis" and "FCgamma receptor phagocytosis", while "regulation of T cell migration" and "cell response to TGFB" were enriched biological processes for mEV miRNA targets.
Discussing more in detail the specific cargo message for colosEVs, we highlighted bta-mir-2284o, among the specific miRNAs, a typical bovine mammary gland miRNA known to be expressed in colostrum and decreased in abundance over time 70 .Also bta-mir-2285cp was exclusively found in colosEVs, which is a liver bovine miRNA that was predicted to regulate amino acids transportation by targeting solute carrier family 7, member 8 (SLC7A8) gene 71 .Another colostrum specific miRNA, reported to play a regulatory role in the immune system, is bta-mir-301a, for which an increased level was observed following T cell activation.This miRNA can activate "nuclear factor kappa-light-chain-enhancer of activated B cells" (NF-κB) signaling, increasing pro-inflammatory cytokines such as interleukin-8 (IL-8), interferon beta (IFN-β), nitric oxide synthase 2A (NOS2) and cytochrome oxidase subunit 2 (COX-2) 72 .Moreover, estrogen receptor α (ERα), one of the estrogen hormone-activated transcription factors, which regulates a large number of genes and is involved in the mammary gland development, is a target of miR-301a 73 .Mouse orthologue of bta-mir-204 can promote the synthesis of milk lipids in mammary epithelial cells by targeting SIRT1 74 .Moreover, in goat, mir-204 together with mir-211 regulates αS1-casein and β-casein synthesis via targeting αS1-casein coding gene in mammary epithelial cells 75 .In bovine, mir-204-5p may target several genes with roles in the nutritional regulation of gene expression in the mammary gland 76  www.nature.com/scientificreports/breast milk promotes the differentiation process of the brown adipose tissue through an expression improvement of the "Runt-related transcription factor 1, translocated to 1" (RUNX1T1) and induces, together with others, the inhibition of the cell proliferation, downregulating the expression of target genes in the p53 pathway 77 .Among specific colostrum miRNA, in a study aimed at investigating the potential regulatory role of miRNAs in the development of gastrointestinal tract, during the early life of dairy calves 78 , mir-211 was predicted to be related to gut epithelial cells and immune cell development, and to inflammatory response.
In the present study, from the 15 functional groups of biological processes enriched for clusters of colosEV miRNA, it has been highlighted processes linked to biosynthesis of macromolecules and regulation of immune signals as "negative regulation of cellular amide metabolic process", "regulation of NIK/NF-kappaB signaling" for cluster 5 or "glycogen metabolic process" for cluster 3, "cellular response to nutrient" for cluster 8 and "macrophage chemotaxis" for cluster 6, as well as several GO terms related to signal transduction in clusters 4 and 5. Interestingly, for GO of upregulated miRNAs in mEV, and thus downregulated in coloEVs, enriched terms related to milk component synthesis were found for clusters 4, 5, 7 and 8, indicating the lack of repression of genes involved in these processes for colosEVs.Moreover, many transcription factors such as SP1 were found as central nodes for both colosEV and mEV miRNAs, indicating a great potential in gene expression regulation by miRNAs once EVs are taken up by recipient cells.Table 3.The restricted list of specific and shared miRNAs in colosEV and mEV cargos used for the functional analysis.
Hub of cluster 1, SNW1, is a protein previously shown to be involved in both splicing and transcription; its role involves binding to the NF-kappaB-p-TEFb complex to facilitate transcriptional elongation of some NF-kappaB target genes.Indeed, SNW1 complex has been identified as a novel transcriptional regulator of the NF-κB pathway 79 .
For cluster 2, several hubs (MED4, MED12, MED30, ARC70, ARC34, ARC77 and ARC205) are genes encoding for components of the mediator complex (MED) also known as activator-recruited cofactor.Mediator complexes are large multiprotein units that regulate gene expression in all eukaryotes and are involved in the transcriptional elongation and termination, mRNA processing, noncoding RNA activation, super enhancer formation, and epigenetic regulation 80 .MED consists of 31 subunits (MED1-MED31) in which MED12 is a critical transducer of regulatory information essential for organogenesis.MED4, instead, encodes for the vitamin D receptor-interacting protein (DRIP) complex, which works as a nuclear receptor coactivator essential for vitamin D metabolism.Moreover, recent studies linked MED4 to epigenetic regulation through a non-canonical pathway where MED4 depletion results in profound changes in the three-dimensional chromatin architecture in contrast to the canonical function of the Mediator complex 81 .Vol:.( 1234567890
In conclusion, the whole message carried by the EVs from the two matrices (colostrum and milk) appears to be similar; the main difference is made by the amount since EVs are 10 to 100-fold higher in colostrum than in milk.
ColosEVs carry molecules, especially miRNAs, potentially capable of modifying metabolic processes of recipient cells involved in signal transduction, cell cycle and immune response, as for EVs of other previously characterized species, but with a special enrichment for miRNAs with epigenetic regulation capacities.
These beneficial characteristics of colosEVs are essential for the calf and could also be exploited for therapeutic purposes in humans, although further studies are necessary to measure the sanitization treatment impact on EV conservation, especially in buffalo where milk is consumed almost exclusively after processing.

Figure 1 .
Figure 1.Morphological characterization of colostrum extracellular vesicles (colosEV) and milk extracellular vesicles (mEVs): (A) Exoview results obtained from the average of the analysis via interferometric imaging (IM) for the CD9 capture spot for colosEVsand mEVs; (B) Transmission electron microscopy (TEM) showing EVs enriched from colosEVsand mEV.Scale bar 500 nm.

Figure 3 .
Figure3.Exploratory analysis on small RNA normalized expression levels for buffalo EVs from colostrum (colosEV) and milk (mEV): (A) Results of the principal component analysis on first 100 features showing that 86% and 6% of the variance is explained by PC1 and PC2, respectively (red dots for colosEV and green dots for mEV samples); (B) heatmap of the seven subjects sequenced for small RNA cargo of colosEV and mEVs with samples in columns and features in rows.

Figure 5 .
Figure 5. Clusters of protein-protein interaction (PPI)-network generated from targets of: (A) colosEV-specific miRNAs (yellow clusters); (B) mEV-specific miRNAs (blue clusters) and (C) shared miRNAs between colosEVs and mEVs (green clusters).Red nodes correspond to central nodes (proteins with more than 50 interactions) and relative protein names are reported close to each cluster (red names indicate miRNA targets).Clusters hereafter are named with progressive numbers as they appear in this figure and with the central node/nodes label.
) was quantified through the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and quality tested by the Agilent 2100 Bioanalyzer RNA assay (Agilent technologies, Santa Clara, CA, USA).The checked RNAs (260/280 and 260/230 ratios and the RNA integrity number) were used to set up a small RNA library suitable for Next Generation Sequencing using QIAseq miRNA library kit (QIAGEN, Germantown, MD, USA) following the manufacturer's instructions.Final libraries were checked with Agilent Bioanalyzer DNA assay and sequenced on a NextSeq 500 (Illumina, San Diego, CA, USA) system, generating paired-end 150 bp sequencies.

Table 2 .
Alignment rate for second mapping step on the Bos taurus genome."C" colostrum EVs; "M" milk EVs.

Table 4 .
Filtered targets based on the number of miRNA hits used for functional analysis.

Table 5 .
First 15 functional groups (lower adjusted p-value) of biological processes enriched for cluster of colosEV miRNA targets.For Clusters 11, 12 and 15, no statistically significant enriched terms were found.

Table 6 .
First 15functional groups (lower adjusted p-value through Benjamini-Hochberg correction) of biological processes enriched for cluster of mEV miRNA targets.For Clusters 7, and 9-15, no statistically significant enriched terms were found.

Table 7 .
First 15functional groups (lower adjusted p-value through Benjamini-Hochberg correction) of biological processes enriched for cluster of shared miRNA targets.For Clusters 4, 5, 12, 13 and 15, no statistically significant enriched terms were found.